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Abstract 

Significant advances, both in the theoretical understanding of rupture processes in heterogeneous 
media and in the methodology for characterizing critical behavior, allows us to reanalyze the evidence for 
criticality and especially log-periodicity in the previously reported chemical anomalies that preceded the 
Kobe earthquake. The ion {Cl~ , K^, Mg^^, and S0^~) concentrations of ground-water issued 

from deep wells located near the epicenter of the 1995 Kobe earthquake are taken as proxies for the 
cumulative damage preceding the earthquake. Using both a parametric and non-parametric analysis, the 
five data sets are compared extensively to synthetic time series. The null-hypothesis that the patterns 
documented on these times series result from noise decorating a simple power law is rejected with a very 
high confidence level. 
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Since the inception of seismology, the search for rehable precursory phenomena of earthquakes has 
shown that the seismic process is preceded by a complex set of physical precursors. In addition to the direct 
seismic foreshocks |2| and seismic precursors ^ , many other anomalous variations of various geophysical 
variables such as electric and magnetic fields and conductivity H, ^ , as well as chemical concentration [0, |8| 
have been documented. However, there is no consensus on the statistical significance of these precursors 
and their reliability ^] , due to a lack of reproducibility and of understanding of the underlying physical 
mechanisms. 

One of the most striking reported seismic precursory phenomena are the time dependence of ion concen- 
trations of ground- water issuing from deep wells located near the epicenter ||^] and the ground water radon 
anomaly preceding the earthquake of magnitude 6.9 near Kobe, Japan, on January 17, 1995. Since the 
first quantitative analysis of this data, which suggested a discrete scale-invariant time-to-failure behavior 
[0], significant advances, both in the theoretical understanding of rupture processes in heterogeneous media 
and in the methodology needed to characterize critical behavior, permits a reassessment of the data. 

Within the critical earthquake model |12, ^, a large earthquake is viewed as the culmination of a 



cooperative organization of the stress and of the damage of the crust in a large area extending over distances 



several times the size of the seismic rupture |^]. Most of the recently developed mechanical models |13, 14, 
H^, IC, 17] and experiments [p!8| , ^] on rupture in strongly heterogeneous media (which is the relevant regime 
for the application to the earth) view rupture as a singular 'critical' point [^]: the cumulative damage D, 
which can be measured by acoustic emissions, by the total number of broken bonds or by the total surface 
of new rupture cracks, exhibits a diverging rate as the critical stress Uc is approached, such that D can be 
written as an 'integrated susceptibility' 

D^A + B{a,-aY, (1) 

The critical exponent < z < 1 is equal to 1/2 in mean field theory |21, and can vary depending on, 
e.g., the coupling corrosion and healing processes. In addition, it has been shown 22, ^ that log- 
periodic corrections decorate the leading power law behavior (||), as a result of intermittent amplification 
processes during the rupture. They have also been suggested for seismic precursors [^]. This log-periodicity 
introduces a hierarchy of characteristic time and length scales with a prefered scaling ratio A [^^. As a 
result, expression (|^) is modified into 

D K A + B{(Tc- af + C {ac- af cos (27r/ In (ctc - a) + (j)) , (2) 

where / = l/ln(A). Empirical |18|, numerical |22,l2^ as well as theoretical analyses [p3 point to a prefered 



value A ~ 2.4 it 0.4, corresponding to a frequency / ~ 1.2 it 0.25 or radial frequency uj = 27r/ ~ 7.5 ± 1.5. A 
value for A close to 2 is suggested on general grounds from a mean field calculation for an Ising and, more 



generally, a Potts model on a hierarchical lattice in the limit of an infinite number of neighbors |24|. It 
also derives from the mechanisms of a cascade of instabilities in competing sub-critical crack growth |22|. 
Empirically, we see that there is not rarely a bias towards a value twice this, corresponding to a better of 
signal-to-noise ratio for rescalings of A^ in the underlying renormalization group equation 

The physical model underlying our analysis posits that the measured ion {Cl~ , K~^, Mg~^^, NO^ and 
SOl ) concentrations of ground- water issued from deep wells located near the epicenter of the 1995 Kobe 
earthquake are proxies for the cumulative damage preceding the earthquake. In this reaction-limited model, 
any fresh rock- water interface created by the increasing damage leads to the dissolution of ions in the carrying 
fluid that can be detected in the wells. We thus expect that the time evolution of measured ion concentration 
should follow closely the equation (|2|). Due to the large heterogeneity of rocks, this 'deterministic' signal 
should be decorated by noise with different realizations for each ion originating from different rock types. 
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To test our hypothesis, we analyze the five ion data sets, thus increasing the statistical significance over our 
previous report However, it is possible that their noise realizations are not completely independent as 
some of the anions and cations are necessarily coupled pairwise. As we shall see, the poor result obtained 
from the analysis of might be due to the fact that it is generally coupled to both Cl~ and NO^ . 

Figure ^ shows the five data sets on which a moving average using three points has been applied. In this 
moving average, the middle point as usual carries double weight except for the two endpoints. Each of the 
five data sets is fitted with equation (|2|) . The time intervals used in the fit were determined consistently for 
all five data sets by identifying the date of the lowest value of the concentration and using that date as the 
first data point. The last data point was taken as the last measurement before the earthquake. This gave 
53 points for the first three data sets, 52 for the fourth and 58 points for the fifth data set. For three of the 
data sets {Cl~ , , and NO^), the fit shown is also the best fit. In the case of SO^ , the best fit had an 
angular log-frequency w « 19 well above the expected range 6 < a; < 9, while the second best fit has to ^ 7 
within the expected range and is shown instead. This is also the case for Mg~^^, where the best fit has a 
very low angular log-frequency w ~ 2 which capture nothing but a slowly varying trend. Here, the angular 
log- frequency of the second best fit a; ~ 16 is approximately double of what is found for the other four data 
sets corresponding to a squaring of A as previously discussed. We note that the ranking of the fits is done 
using the variance between the data and the fit with equation (|2|) , where the stress a has been replaced with 
time t, tc being the critical date of the earthquake. This assumes a Gaussian white noise-distribution, which 
constitutes a reasonable null-hypothesis for the noise but presumably does not fully refiect the reality. 

We observe that the values of the critical time tc predicted for the earthquake from the various fits are 
not only very stable but also remarkably consistent with the true date 95.053 of the earthquake. And except 
for the case of Mg~^^ , where it has approximately doubled, the same is true for the angular log- frequency lv 
(although in 2 cases, we took the second best fit as explained above). This robustness with respect to the 
date of the earthquake and the preferred scale ratio A = e'^'^^^ ~ 2.5 is quite remarkable, considering that 
the value obtained for the exponent z varies by approximately a factor of 4 from the smallest to the largest 
value. 

We complement these fits by a direct 'non-parametric' analysis of the log-periodic component, taken as a 
crucial test of the critical earthquake model captured by (|l]) , (|2|) . After de-trending each of the five chemical 
time series c{t) using 

^ ' c{tc-ty ' ^ ^ 

which should leave a pure log-periodic cosine if no other effects were present, we apply a Lomb periodogram 
p^ ] to the de-trended data as a function of log {{tc — t) /tc)- Lower fight figure of |l| shows the five spectra 
superimposed. In all case except for , we observe significant peaks with a log-frequency / = ^ between 1 
and 2, i.e. u; is approximately between 6 and 12.5, and in two cases {Cl~ and NO^), we have two very clear 
peaks at / ~ 1.2 {to ~ 7.5) standing out with a Lomb weight of ~ 11 and ~ 13, respectively, in agreement 
with our prediction. 

In order to assess the significance of these results, we present rather exhaustive statistical tests performed 
by constructing synthetic time series that differ from the real data only by the log-periodic pattern and we 
of course follow the same testing procedure. 

There are several reasons why the results from the analysis above should be compared with those of 
synthetic tests. In particular, the real time series have been sampled non-uniformly in time because the 
water from the wells was collected from commercial bottles of which the production dates could be identified. 
We observe that the depletion process of bottles in stores implies that the number of bottles with a production 
date t prior to the date of the earthquake tc is proportional to log (tc — t), leading to a uniform sampling 
in log {tc — t) instead of time tc — * to the earthquake. In order to investigate the effect of this sampling of 
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Figure 1: The concentrations of CI (top left), (top right), Mg^'^ (middle left), A^'Og (middle right) 
and S0^~ (lower left) as a function of date. The fit is equation (|2|). The values of the fits are Cl~ : 
A « 15.3, B -1.9, C « 0.34, z « 0.58, tc ~ 95.060, ^ -0.40 and w 7.0. K+: A ^ 0.31, B ^ -.044, 
C ^ 0.036, z ^ 1.2, tc ~ 95.061, (/> ^ 0.05 and w 7.2. A 5.5, B ^ -0.65, C 0.11, z « 0.82, 

« 95.069, ?a -1.3 and lu ^ 16. iVOg": A ^ 11.6, B « -3.7, C ^ -0.67, z ^ 0.28, tc ~ 95.059, 2.1 
and w w 7.5. SOJ": ^ 17.3, B ^ -2.2, C ^ 0.83, 2 0.77, tc ~ 95.066, « -.31 and uj « 7.4. 
Lower right figure is the Lomb periodogram of the five time series as a function of log-frequency defined as 
the conjugate to log ((tc — t) /tc). Prior to the spectral analysis, the leading power law has been removed 
using equation (P). The values of the Lomb peaks are 11, 3, 6, 13 and 5 in the same order as above. 



the Kobe data with respect to log-periodic signatures, we have as a first step generated twenty synthetic 
data sets each with 56 points as in the original data analyzed in ref.|ll]. The synthetic data was generated 
using a noisy power law with the parameter values of the leading power law of the fit presented in ref.|ll|. 
Specifically, the equation c (U) = 15.45 + (—1.877 + k ■ (0.5 — ran)) • (95.053 — U)^'^'^ was used with different 
noise-levels k. Here, ran is a uniform random number generator with values in the interval [0 : 1[ from |26|. 
In order to obtain a sampling similar to the one found in the original data, the sampling times tj of the 
synthetic data was chosen as tj = 95.053 — (95.053 — 93.847)e~^'"^'^, where 93.847 was the date of the first 
point in the original analysis. This then gives us a distribution of sampling times, which are uniform in 
log-time to t^, where t^. = 95.053 is the time at which the Kobe earthquake occurred. As noise-level, we 
used k = 0, k = 0.37 and k = 0.73, corresponding to no-noise, a noise level of the same amplitude as the 
estimated log-periodic oscillations and a noise level twice that. 

In the case of no noise k = 0, the Lomb Periodograms exhibit in general a forest of peaks over a large 
log-frequency interval which average out to give a flat spectrum. Adding noise of amplitude k = 0.37 
and k = 0.73 enhances the peaks but the averaging over the periodograms of the different synthetic data 
sets again produces an essentially flat spectrum. By "flat", we mean that, in each individual synthetic 
time series, about 7 "peaks" are found above 70% of the largest one over the entire log-frequency interval 
< / < 6. This suggests a qualifier for the significance of a Lomb peak, corresponding to counting the 
number of crossings between the periodogram and a 70%-level or a 50%-level, taking the highest peak as the 
100% reference (a single peak corresponds to two crossings). For all cases, we obtain at least 10 crossings 
depending on how one defines a peak. We also see that the noise has the effect that the largest peak is no 
longer to be found for the lowest frequencies. 

We now prodeed with synthetic tests. As a first null-hypothesis, we take white noisy power law. An 
indication of the statistical significance of the results obtained from the chemical data can be estimated 
by comparing the periodogram of the de-trended data (|3|) with the periodograms of a sequence of random 
numbers uniformly distributed between and 1 after a moving average over three points has been applied. 
Hence, we have generated 1000 sequences with approximately the same number of points (53) and the time 
interval of the true data: none had a Lomb peak above 10. Out of 10.000 synthetic data sets, 9 had 
Lomb peaks above 10, but none above 12. An additional 100.000 data sets were generated providing three 
Lomb peaks above 12, but none of them had a log-frequency between 1 and 2. For peak values above 
10, the 100.000 data set had six such peaks in that log-frequency range. If we use a Gaussian white noise 
distribution instead of the uniform one, out of 1000 synthetic data sets, we get 20 peaks with a log-frequency 
between 1 and 2 and a Lomb weight above 6, but none above 10 with a log-frequency between 1 and 2. Out 
of 10.000 data sets, only one such peak is observed. 

Many other systematic and more elaborate tests have been performed with and without smoothing, as 
well as by directly generating noisy power laws and then performing the fit with (^) and the de-trending 
with (^). We find that smoothing after de-trending is not significantly different from smoothing before 
de-trending. Similarly, a generation of noisy power laws gives results similar to the foregoing white noise 
hypothesis. 

Figure ^ presents a global summary of our tests by showing the bivariate distribution of angular log- 
frequencies uj and Lomb peak h obtained from a series of 100,000 synthetic time series. The five vertical 
lines correspond to the five real times series and can be obtained from the lower right panel in figure 1 by 
locating the log-frequency and Lomb peak height of the highest peak for each data set. One can visualize 
that not only the fundamental frequency but its harmonics are sometimes visible {Cl~ and NO^) and, as 
we said above, may be dominating (Mg~^^ and S0^~). 

The interpretation of these results is that there is a confidence level of 99.99 % for a single peak with 
a log-frequency between 1 and 2 and a Lomb weight above 10. This confidence interval is evaluated with 
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Figure 2: Bivariate distribution for the null hypothesis of angular log-frequencies u and Lomb peak h 
obtained from a series of 100,000 synthetic time series. The five vertical lines correspond to the five real 
times series. They are obtained from the lower right panel in figure 1 by locating the log-frequency and 
Lomb peak height of the highest peak for each data set. 
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respect to our initial null hypothesis of of uncorrelated (white) noise and would of course change with a 
different null hypothesis. Nevertheless, the confidence level achieved with the present null-hypothesis is 
above 99.9% and we hence expect it to remain would remain very high over a range of null- hypothesis. Also, 
the confidence level will also change depending on whether one assumes that the various ion concentrations 
are independent or not. 

In conclusion, we wish to stress that the presented analysis constitute but a single case-study. Hence, it 
does not propose a recipe for earthquake prediction. However, we feel that the statistical evidence for this 
particular analysis is significant enough to encourage further studies along similar lines. 

Acknowledgments: We are grateful to H. Wakita for kindly providing the data and for useful correspon- 
dence. H.S. thanks Y. Huang for collaboration at an early stage of this work and for many useful discussions. 
For further discussion of statistical tests and interpretations from a different point of view, see |27]. 
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